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Abstract 

We study the consistency of the description of charge distributions and radii of nuclear clusters 
obtained with semiclassical nuclear pasta models. These nuclei are expected to exist in the low 
density outer crust of neutron stars. Properties of the arising clusterized nucleon matter can be 
compared to realistic nuclear properties as experimentally extracted on earth. We focus on non 
iso-symmetric light clusters with nucleon number 8 < A < 30 and use Monte Carlo many-body 
techniques. We simulate isotopic chains for a set of selected nuclei using a model Hamiltonian 
consisting of the usual kinetic term, hadronic nucleon nucleon (NN), Coulomb and an effective 
density dependent Pauli potential. It is shown that for neutron rich (deficient) clusters neutron 
(proton) skins develop. Different (matter, neutron, proton, electric charge) radii are computed for 
this set of non iso-symmetric nuclei. Nuclear binding energies are also analyzed in the isotopic 
chains. 
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I. INTRODUCTION 



The knowledge of nucleon distribution in finite nuclear systems is a relevant topic in nu- 
clear physics since it is not only important in characterizing the spatial extent of the nucleus 
but it is related to other issues like, for example, the location of the nucleon drip-lines |l| 
or the determination of nuclear masses in the r-process nucleosynthesis in the astrophysical 
context . A variety of approaches have provided with nuclear mass formulae that allow 
testing properties of currently available nuclei from experiment Both non relativistic 
density-dependent Skyrme js] and Gogny interactions ^\ and relativistic mean field (RMF) 
calculations [tI have been largely used as they provide a successful reproduction of charge 
distribution and size of nuclei. Other many-body techniques like Monte Carlo, classical 
Molecular Dynamics, Antisymmetrized Molecular Dynamics or even more elaborate as the 
auxiliary field diffusion Monte Carlo have the potentiality of being applied to this type of 

However, in the context of direct simulation of cluster- 



calculations 
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ized matter in systems such as the outer crust of neutron stars or supernova matter these 
techniques are not suitable. Due to the nature of these intensive calculations it is usual to 
construct nuclear interaction models in a way so that they will retain the essential physics 
involved, yet they keep as simple as possible for saving CPU time. 



In exploratory works 
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13l | we used Monte Carlo techniques to study nuclear binding 
energies and nuclear density distributions for iso-symmetric spin saturated medium mass 
nuclei with nucleon number up to A 100 using semi classical simulations via Monte Carlo 
many-body techniques. In [ij] the purpose was to study the role of an effective Pauli po- 
tential dependent on density parameterized by one nuclear variable, the Fermi momentum 
of the constituent nucleons for nuclei with A up to ~ 50. In this semiclassical simulation 
of a finite nucleon system it was demonstrated that the empirical binding energies for these 
nuclei can be reproduced satisfactorily even if one uses a simplified NN interaction provided 
a density dependent Pauli potential is included. In that case a delicate counterbalance of 
kinetic and potential energy contributions arises allowing for the nuclear saturation mech- 
anism. In particular, only isospin and additionally spin dependence in the NN interaction 
sector was included to size the impact on binding of (anti) alignement of nucleon spins. Ob- 
tained results show that this latter improvement in the description of the interaction weakly 
impacts the reproduction of experimental nuclear binding energies. 
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In [l3| we considered a wider range of nuclei with 8 < A < 100 for spin-isospin saturated 
Z = N (even Z and N) nuclei. We showed that as the nucleon Fermi momentum is no 
longer a good parameter due to the fact that the average Fermi momenta in nuclei with 
A > 50 saturate to a value ~ 260 MeV/c (e.g., 265 MeV/c for a ^'''^Pb nucleus) Q, and 
additional dependence must be considered in order to parametrize fermionic correlations. 
This result, then, generalizes the previous conclusion that the density dependence of the 
Pauli potential is crucial to reproduce the empirical nuclear binding energies [l2 1 . If heavier 
systems are explored, one obtains that for systems with nucleon number A larger than ^ 200 
the binding energies per particle and Pauli potential strength contribution stabilizes and the 
corresponding binding energy per particle for symmetric nuclear matter at low temperatures 
of ~ —16 MeV is obtained. 

As mentioned, in this type of semiclassical simulations one considers effective degrees 
of freedom under the form of nucleons. In this sense this treatment is appropriate to the 
interaction range of fm (or energies of the order of the pion mass m^r ^ 140 MeV) in 
an analogous way to the very successful effective field theories (EFT) used in hadronic 
physics |l5l ]. 

Other type of nuclear systems of interest in fields such as nuclear astrophysics can 
be considered in this simulation scenario. Clustered configurations arising in low den- 
sity neutron star matter, below densities ~ O.lpo (Po = 0.17 fm^^) and known as pasta 



phase 



Q, Q, Q, 



19| can be studied using this approach, beyond the usually adopted 



Wigner-Seitz cell approximation that lacks the inclusion of higher order correlations 20|. 
In simulations of heavy-ion collisions the fragmentation yields at late times show that the 
final mass distribution of clusters is a relevant topic that largely depends on higher order 
correlations in the simulations [21]. 

From the experimental point of view there is an increasing effort in the experimentation 
with radioactive beams such as, among others, ISOLDE I22I, those at RIKEN 23] or the 
proposed RIA [2J] and FAIR [25| . Special attention is devoted to the new available (proton) 
neutron rich nuclei where observables such as the (proton) neutron skin radius can be mea- 
sured. Experiments, like e. g. P-ReX at Hall A in Jlab, will measure the neutron radius 
in 208 using parity violating electron scattering j26|. Additional phenomenology such as 



halo nuclei 



271 ] has also been observed. 



The structure of the manuscript is as follows. In section [TTl we present the effective nuclear 
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model used in our simulations. In section IIIII we present the results for the observables 
obtained from the calculations of the different nuclear isotopic chains considered in this 
work and finally we conclude in section [IVl 



II. NUCLEAR SIMULATION MODEL 

In the present approach we describe matter in nuclear astrophysical scenarios treating 
nucleons as classical, structureless particles. The model Hamiltonian consists of nucleon 
kinetic energy, nuclear interaction (Vatat), Coulomb (Vcoui) and Pauli {Vpaun) potentials. 
The Pauli potential depends not only on relative positions but on momenta of the interacting 
nucleons. It simulates nucleon fermionic nature using a Gaussian functional form introduced 
by Dorso et al. [28] allowing, however, for density dependence, see below. 

The model Hamiltonian is given by, 

A p2 A 

H = J27r^+ J2 ["^NNirij) + Vcoui{rij) + Vpauii{rij,Pij)], (1) 

where pi is the 3-momentum of i-th nucleon and Vij = |ri — r^l {pij = |pj — Pj|) the 
relative distance (momentum) of the i-th and j-th nucleons. m-Ar is the nucleon mass. We set 
h = c=l. 

Since we are interested in the description of finite clusters we will use a simplified phe - 
nomenological NN potential interaction, similar to those used in nuclear pasta models 19| . 
consisting of two Gaussians as follows, 

VNNin,) = + [b + cr,r,]e-^-/"^ (2) 

where Tj {tj) is the isospin third-component of i-th (j-th) nucleon (-|-1 for protons, —1 for 
neutrons). We will label this interaction potential as GM model. Let us emphasize that this 
model has a functional dependence that, although simple, can be used to model extended 
nuclear systems at high densities (several times nuclear saturation) 1^. It includes an 
isospin dependent part to ensure that while pure neutron matter is unbound, symmetric 
nuclear matter is appropriately bound However, in this work we are mainly interested 
in studying the low density limit of nuclear systems retaining the same functional form in the 
NN interaction. In this way a consistent check can be performed using the same treatment 
for nuclear systems in the whole density range where the nucleon picture holds. It has been 



shown, for example, that the lack of consistency between results for matter response to a 
weakly interacting probe, as obtained with fine tuned models for a limited range of systems, 
can lead to, for example, overpredict neutrino opacities in the astrophysical context 19|]. For 
accurate nuclear shell calculations one should relay, however, on more elaborate treatments 
including the before mentioned richer NN models with additional two and three-body forces 
as in Quantum Monte Carlo 29| or ab initio calculations [30|. 



TABLE I: Model parameters. 



a{MeV) b{MeV) c{MeV) Ai (/m) As (/m) 



50 



0.8 



In order to calibrate the model we choose a set of mirror nuclei where the neutron number, 
A^, and the proton number, Z, can be interchanged as, for example, in ^^C and ^^O and fit the 
experimentally available root mean square (RMS) matter and charge radius. The parameter 
values used in this work are given in Table [11 Proton and charge radius can currently be 
measured with high relative precision using electromagnetic probes by the cross section 



measurement of elastic scattering of electrons 311] or muonic atoms 



32]. However, for matter 



and neutron radii the accuracy is not so good due to the fact that there is uncertainty in the 
nuclear regions probed by the various processes and due to the effective nuclear interactions 



331 ] provide information on neutron distribution in 



between nucleons. Antiprotonic atoms 
a nucleus. 

Since protons are electrically charged particles we introduce the Coulomb potential, 

Vcou^in,) = (3) 

where e is the proton electric charge. 

The antisymmetrization is obtained in an effective way by using a Pauli potential under 
the Gaussian form given by, 

2 2 
r . p . 
»J 'J 

Vpauuir-ij.pij) = Vp e (4) 

where dr^Tj {^a^aj) is the Kronecker's delta for the isospin (spin) third-component. It prevents 
nucleons from occupying the same phase space volume when they have the same quantum 
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numbers. We will allow a density dependence for this Pauli potential in the parameters 
qo,Po and Vp as in [12I, Although this effective treatment using a Pauli potential is a 
way of providing antisymmetrization it has some well known shortcomings 34 1. 

In our simulation a nucleus is considered to be a composite object where A non-relativistic 
point-like nucleons aggregate according to the interaction given by the Hamiltonian in Eq.([T]). 
The typical internucleon distance in an iso-symmetric nucleus is given by ~ 2r = ( 
where p = 2p|./37r^ is the nucleon density and the nucleon Fermi momentum. In Moniz et 
al. Fermi momentum is asignated by using a simplified Fermi gas model to parameterize 
the density dependence of electron scattering cross section on a set of selected nuclei. To 
study the structure of the clustered configurations arising in the simulation one can compute 
the RMS radius by averaging the corresponding density operator over the simulation box 
volume (y = L^). 

In this work we will consider different charge (matter, neutron, proton) density operators 
that will be designated as pi with i = {m, n,p} and charge qi with i = {1,1/2 — Tj, 1/2 + r,} 



respectively and j 



A. The charge density operators can be written generically as 



(5) 



Let us note that the RMS radius is computed by averaging the corresponding charge operator 
over the box volume. In this way the RMS radius can be obtained as 



< rt > 



1/2. 



/ d^rpi 



r r 



1/2 



(6) 



Jd^rpi{r) 

in order to compare with the nuclear charge radius, Vch, that can be experimentally measured. 



it is usual 



35| to assume a Gaussian charge distribution inside the proton so that an effective 



proton finite size, Sp, is used and typically Sp 
for the charge RMS radius is. 



0.8 fm 



The expresion that will be taken 



l/2_ 



211/2 



(7) 



where < >^/^ is the point proton charge RMS radius. 

Having explicited the different RMS radii above we can proceed for neutron rich (deficient) 
clusters and define the difference of RMS radii of the neutron and proton distribution js^] , 

^^npi as 

(8) 



An 



np 



^ I p ^ 



In a similar way we define Arp„ = — Ar„p for tfie proton ricli clusters. There are estimates 



of this observable in the literature such as in 



38| where they quote the droplet model 



based value ARMS = |t where t is the neutron skin thickness or the shell model potential 



extimates 39|]. Note that diffuseness of the surface and non-uniformity of densities add 
corrections to this last expression. This observable is also closely related 401] to the shift 
in the surface to volume ratio of symmetry energy, rs/v = O's/O'V as derived in Skyrme- 



Hartree-Fock models 



4l[ | and other works 42|, |43 1 . 



In [4J] a linear fit is given using a set of experimental data measured in antiprotonic atoms 
as the isospin asymmetry parameter / = [N — Z)/A variates. They give an expression for 
the fit of the neutron skin 



Ar„p = (0.9 ± 0.15)/ + (-0.03 ± 0.02) fm 



(9) 



Notice however that this isjust a fit to medium size nuclei. For smaller A there are also 



some calculations such as 



45 



461] where they consider the light nuclear sector. Additional 



information using this data can be derived for the symmetry energy [47]. 



III. RESULTS 



In the simulation we initially place the finite nuclear system under study, consisting of 
A nucleons distributed uniformly inside a sphere of radius Rq in the range 2 — 3 fm, within 
a cubic box of volume V = Li^ and impose L » Rq. Computer simulations are performed 
by Monte Carlo techniques where using the Metropolis algorithm 48| the ground state 
configuration is searched by thermal relaxation at a fixed temperature T = 1 MeV. Note 
that the role played by temperature in this semiclassical approach is just that of a parameter 
to somewhat simulate the effect of zero-point motion in a phenomenological way. Once the 
thermalized stable clustered configuration is obtained, a nucleus, we statistically sample the 
system for the averaged magnitudes, namely, binding energies and RMS radii. ^ ^ 

In Fig. [U we can see the energetic contributions to the Mg isotope chain [49] . Solid, 
dotted and dashed lines show the potential, kinetic and total binding energy, respectively. 
Positive (kinetic) and negative (potential) energetic contributions allow to reproduce the 
experimental isotopic binding energy values. Note that due to the Pauli potential the kinetic 
energy per particle departs from the value k = ^ksT as one would expect in an uncorrelated 
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classical plasma. 
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FIG. 1: Potential (solid line), kinetic (dotted line) and total (dashed line) binding energy per 
nucleon for the isotopic chain of Mg. 



In Fig. [2] we show the shift in binding energy per nucleon {6E) from the spin saturated 
iso-symmetric species (A^ = Z) versus nucleon number A for Carbon (C), Oxygen (O), Neon 
(Ne) and Magnesium (Mg) isotopes as obtained by our simulation. We can see that for Ne 
and Mg the minimum is shifted to the neutron rich side in agreement with experimental 
binding energy data j^]. 
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FIG. 2: Energy shift per nucleon with respecto to the N = Z species for the isotopic chains of C, 
O, Ne and Mg nuclei. 



In Fig. [3] we plot matter (boxes), neutron (circles) and charge (triangles) RMS radii 
for the Carbon isotopic chain versus the nucleon number A. We can see that for the iso- 
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symmetric species with Z = 6, N = 6 there is no skin present. As A increases this neutron 
skin grows since the charge radius remains approximately constant. For the proton rich 
isotopes a proton skin develops while for nuclei on the neutron rich side a neutron skin 
forms in agreement with more elaborated relativistic mean field jsc)] or relativistic Hartree- 
Bogoliuvov [46] descriptions including nuclear deformations. 
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FIG. 3: From top to bottom: neutron, matter and charge RMS radii for the isotopic chain of 
Carbon. 



We can see more in detail in Fig. H] matter (boxes) RMS radii for the Carbon isotopes 



along with experimentally deduced values from Ref . 5l(| • Data are extracted from hadronic 



reactions involving break up of nuclei. For example, in the reaction 

a + A^b + X (10) 

the width of the momentum distribution, a, depends on the mass number of the projectile 
(Aa) and of the fragment (Ah) but it is found to be essentially independent of the target 
mass and beam energy. It can be expressed using the Goldhaber model |5^] as 



a 



^jA^{Aa - A^)l{Aa - l)ao (11) 



with this prescription and in the approximation where each isotopic component is considered 
as a Fermi gas with Fermi momentum pp, then the nucleon density is p = and ctq is the 
reduced width of the momentum distribution and it is associated with the Fermi momentum 
of the nucleons, uo = pf/\/^- The interaction cross section, o"/, can be experimentally 
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measured with 5% accuracy and functionally can be parameterized as 



(12) 



Ri is the interaction radius for the selected projectile on the target nucleus. Fitting cor- 
respondingly To to the experimental cross section data one obtains the general trend for 
the interaction cross section. For neutron rich isotopes there is, however, a deviation from 
this prediction as a neutron skin forms. Using Glauber model calculations the matter RMS 
radius can be estimated by using the approximate relation above [27 1. 
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FIG. 4: Matter RMS radii for the isotopic chain of Carbon. Empirical data are shown along with 



error bars taken from Ref. 



5l| 



In Figs. [5] and [6] we plot charge and neutron RMS radii, respectively, for the Carbon 
isotopes along with experimental data. We can see that while the neutron RMS radius is a 
rising function of A the nuclear charge radius stabilizes from the iso-symmetric isotope (T^ = 
0) up to values = 4 where = (A^ — Z)/2 is the third-component of the isospin. Note 
that in our analysis we do not consider deformated nuclei since a many-body r-dependent 
potential used in nuclear pasta calculations is used. Deformations, /32, in the carbon isotopic 
chain range from [—0.4,0.4] 46 1. 

In Fig. [7] we show matter RMS (boxes), neutron (circles) and charge (triangles) RMS radii 
for the Oxygen isotopes. We can see that as isospin asymmetry grows a neutron (proton) 
skin develops on the neutron (proton) rich side, that is thinner than in the case of Carbon 
as shown in Fig. [31 
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FIG. 5: Same as Fig. IHfor the charge RMS radii. 
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FIG. 6: Same as Fig. [Hfor the neutron RMS radii. 



In Fig. [8] we can see the difference in neutron and proton RMS radii, A„p, as defined in 
Eq.Q, versus nucleon number A for the Neon isotope chain, as obtained in more sofisticated 
calculations as those of the RMF approximation with the constant gap approximation in 
the axially deformed oscillator basis (DEF) model from Ref. 53| and compared to our 
semiclassical simulation with the GM model. We can see that the divergence of both models 
is bigger as the isospin asymmetry grows. However if we consider the result from the linear 
fit in Eq. ([9]) at the highest isospin asymmetry considered in this chain, at A=30, the 
discrepancy is about 16%. This shows that the description provided by a simplified nuclear 
interaction, as that of the GM model, similar to those used in nuclear pasta models allows 
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FIG. 7: Matter (boxes), neutron (circles) and charge (triangles) RMS radii for Oxygen isotopes. 



a reasonable consistent description as compared to more elaborate treatments. 
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FIG. 8: Difference in proton and neturon RMS radii versus nucleon number A for the Neon isotope 
chain as deduced from the DEF model 53] and compared to our calculation with the GM model. 



In Fig. [9] we show matter and charge RMS radii versus neutron excess for the Mg isotopic 
chain. Experimental values for matter radii are depicted for comparison as taken from 



Ref. 



511 ]. Neutron excess (fault) in the isotope chain results in a developing neutron (proton) 



skin thinner than in the case of Carbon and Oxygen isotopes. As we can see there is a 
competing effect that tends to enlarge the neutron skin as the isospin asymmetry grows and 
that tends to decrease the neutron skin as the atomic number incresases. 



12 



5 

4.5 

u) 4 

2 3.5 
c/) 

QC 3 
2.5 
2 

-4 -2 2 4 6 
N-Z 

FIG. 9: Matter and charge RMS radii versus the neutron excess for the Mg isotopic chain. Ex- 
perimental values for matter radii are depicted for comparison as taken from Ref. [51,]. 

IV. CONCLUSIONS 

In this work we anlyze the consistency of the description of finite nuclear systems, clusters, 
in the limit of low density of interest in the astrophysically relevant nuclear pasta models. 
In this regime densities in the range of tenths to several times nuclear saturation density 
are considered. We obtain some observable properties of finite nuclei regarding charge 
distribution and different (matter, neutron, proton, charge) RMS radii for spin saturated 
nuclei in the mass range with 8 < A < 30. We use many-body semiclassical simulations 
with Monte Carlo techniques to compare the results provided by such a nuclear interaction, 
designed to be used in many-body systems to experimental data or more elaborated nuclear 
models. 

We use a hamiltonian model where apart from the kinetic term, the NN interaction 
consists of a hadronic. Coulomb and Pauli effective potential terms. We have shown that 
for a set of nuclear isotopic chains (C, O, Ne, Mg) experimental binding energies can be 
well reproduced if a density dependence for the Pauli potential is allowed. We have also 
analyzed other nuclear structure properties such as matter, neutron and charge radii. We 
have found, besides, that our obtained radii agree reasonably well with those deduced from 
experiment. As the isospin asymmetry parameter grows, the departure from this trend, 
as obtained with our simplified model, is larger remaining below 16% however. 

For non iso-symmetric clusters we find that a neutron (proton) excess results in a neutron 
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(proton) skin developing in the exterior shell of the nuclei and computed as the difference of 
neutron and proton RMS radii. For iso-symmetric nuclei in our model there is a vanishing 
skin, in accordance with the liquid drop formula. Charge radii are computed in our simula- 
tion from point-like particles by folding the proton RMS radius with the proton finite size 
Sp ~ 0.8 fm. Our computed values agree with measured values within experimental error 
bars. Charge radii stabilizes as A increases giving the possibility of neutron skin formation 



on the neutron rich sector. However, in particular, the Neon isotope anomaly [53|] can not 
be reproduced within our simplified spherical configuration. For this type of phenomenology 
nuclear deformations play an important role. Additional enrichment of the model using two 
and three-body terms could provide better agreement. Using these techniques of semiclas- 
sical simulations on effective degrees of freedom allows a consistent procedure to further 
testing asymmetric systems by studying the formation of clusters in these environments and 
how their structure affects the opacity of this type of matter. 
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